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A multi-chain mean-field theory is developed and applied to a two-dimensional system of weakly 
coupled S = 1/2 Heisenberg chains. The environment of a chain Co is modeled by a number of 
neighbor chains Cs, S = ±1, . . . ± n, with the edge chains C±„ coupled to a staggered field. Using 
a quantum Monte Carlo method, the effective (2n -I- l)-chain Hamiltonian is solved self-consistently 
for n up to 4. The results are compared with simulation results for the original Hamiltonian on large 
rectangular lattices. Both methods show that the staggered magnetization M for small interchain 
couplings a behaves as M ~ ^/a enhanced by a multiplicative logarithmic correction. 

PACS numbers; 75.10.Jm, 75.40.Cx, 75.40.Mg 



Quasi one-dimensional (ID) quantum spin systems 
have become an important field of study in solid state 
physics. Many unusual, theoretically predicted proper- 
ties of ID systems have been observed in real materials. 
For example, the gapless two-spinon spectrum Q| of the 
S = 1/2 Heisenberg chain has been observed in neutron 
scattering experiments on KCuFa ||^, and the Haldane 
gap predicted for integer S ^ has been detected, e.g., in 
the S* = 1 compound CsNiCls Q . Quantum critical scal- 
ing Q has been observed in the NMR relaxation rates of 
the S ^ 1/2 system Sr2Cu03 possibly even including 
anticipated logarithmic corrections ||^. In spite of 
the success of strictly ID models for these and many other 
quasi-lD magnetic materials, interchain couplings can be 
important as well. A single isotropic chain cannot order, 
not even at T = 0, whereas a transition to a Neel ordered 
state is often observed at low temperature; KCuFs and 
Sr2Cu03 both order at Tjq « 5 K. Interchain couplings 
also change qualitatively the nature of the low-lying ex- 
citations and lead to interesting dimensional cross-over 
phenomena. 

One way to take into account interchain couplings J± 
in a quasi-lD system with long-range order is to model 
the environment of a single chain Cq by a staggered mag- 
netic field §-11. The effective ID system can be solved 
numerically on small lattices [|l^,^, or using analyti- 
cal techniques [|ll[. In this Letter, the mean-field ap- 
proach is extended to include also a number of neighbor- 
ing chains Cs to which Co is coupled. In two dimensions 
S — ±1, . . . ± n. A staggered field is coupled to the edge 
chains C±„, to model their long-range ordered environ- 
ment. Fluctuations neglected in the environment of Cs 
are approximated by a modification of their intrachain in- 
teractions, in such a way that self-consistency is achieved 
in the induced staggered magnetizations on Cq and Cs- 
The effective {2n + l)-chain Hamiltonian can be solved 
using numerical methods, which typically perform much 
better for a few coupled chains than for 2D or 3D lattices. 

Here these ideas will be applied to a system of antifcr- 
romagnetic Heisenberg chains, with the Hamiltonian 



H — J Sij ■ Si+i J- -|- J± Sij ■ S 



(1) 



where Sij denotes a spin-1/2 operator at site i of chain 
j. The focus will be on the dependence of the T — 
staggered magnetization M — (S*! j) (— 1)'^-' on the cou- 
pling constant ratio a = J±/J. The question of whether 
or not long-range order (M > 0) develops for arbitrarily 
small a > has been the subject of numerous studies. 
Conventional spin-wave theory predicts a finite critical 
value ac below which M = but RPA Q gives 

CKc = 0. Some self-consistent calculations predict ac — 
p^ , whereas others have given ac as high as 0.2 [ p^ . 
Renormalization group analyses of the interchain inter- 
actions are associated with subtleties [^7|jl^], and com- 
pletely conclusive results have not been presented; how- 
ever ac = appears most plausible | . An analytical 
treatment of the single-chain mean field theory gave the 
behavior AI ~ -^/a for small a Numerically, M has 
been calculated using exact diagonalization (l^Jl^ and 
series expansion techniques the former indicating 
ac ~ 0.1 — 0.2, and the latter giving an upper bound 
ac ^ 0.02. Numerical calculations have in general been 
hampered by convergence problems and difficult extrap- 
olations for small a. Here multi-chain mean-field cal- 
culations will be complemented by large-scale quantum 
Monte Carlo simulations of the original 2D Hamiltonian 
(|l|) . It will be shown that quadratic lattices are not suit- 
able for extrapolations to the thermodynamic limit when 
a ^ 1, due to unusual, non-monotonic finite-size effects. 
Using rectangular lattices with aspect ratios L^/Ly as 
large as 16 it was, however, possible to study systems 
with a as low as 0.02. Both the mean-field calcula- 
tions and the 2D simulations indicate that M vanishes 
as a ^ slower than ^/a, due to a logarithmic correc- 
tion to this form. 

In the conventional single-chain mean-field treatment 
of the Hamiltonian (0) p-|lT|, the coupling of a chain 
j to its nearest-neighbors j ± 1 is approximated by 
J±J:^Sl,[{Sl^^,) + (Sl^+i)]. In a Neel state (Sl^) = 
{—iy~^^M, and one obtains an effective ID Hamiltonian, 



1 



H, = jY,s,-s,+,-hY,i-^ys!, (2) 

i i 

with the self-consistency condition h — 2J±M which di- 
rectly relates M{h) to M(J_l) of the 2D system. 

The idea of the multi-chain mean-field theory is to 
model the environment of a chain Cq by its first few 
neighbor chains Cs, 5 = ±1, . . . ±n, with only the edge 
chains C±„ coupled to a staggered field. This induces a 
staggered magnetization in all chains. The dynamic en- 
vironment for Co provided by the C's chains should be 
considerably more realistic than just the static staggered 
field of the single-chain theory. If Co and Cs are iden- 
tical chains, it is not possible to obtain a self-consistent 
description, however. The staggered magnetization will 
be largest at the edges and decrease towards the cen- 
ter, due to the neglected quantum fluctuations at the 
edges. These fluctuations can be approximated by a mod- 
ification of the intrachain interactions of Cs- There are 
clearly many possible ways of doing this, and the op- 
timum way, that would best mimic the presence of an 
infinite half-plane of other chains, is not obvious. One 
requirement is that the additional interactions have to 
be invariant under spin rotations in the xy-plane (since 
the field breaks the 0(3) symmetry, an 0(2) symmet- 
ric effective interaction in Cs is permissible). Here the 
simplest interaction satisfying this requirement will be 
considered, namely, the xy part of the coupling is given 
a strength J^^^ ~ J(l -|- \\s\) different from J|^| = J. 
Increasing A|5| > increases the quantum fluctuations. 
The (2n -I- l)-chain effective Hamiltonian is then 

L n L n—1 

Hn ^ X! ' ^^+1 J + X! X! ' 

2—1 j—~n i—1 j——n 

L n 

+ X! X! -^I-^I i^lsSi+i.s + SlsS^i+i^s) (3) 

i=l \S\ = l 
L 

i=l 

There are n -I- 1 self-consistency conditions: 

M = A/o = Afi = M„, 

h = MJ^. (4) 

Since the environment of a chain Ck becomes more sim- 
ilar to that of the real 2D system the closer it is to 
the center [k = 0) of the effective 2n + 1 chain system, 
the self-consistent parameters can be expected to satisfy 
< Ai < . . . < A„. For a given a, the magnetization (as 
well as other properties) should converge to its correct 
value as n,L —^ oo. Therefore, the details of the intra- 
chain interactions used to achieve self-consistency can be 
seen to be unimportant; they will only affect the rate of 
convergence with increasing n. 
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FIG. 1. Quantum Monte Carlo results for the staggered 
structure factor on rectangular lattices with different aspect 
ratios and a = 0.05. The behavior for a quadratic lattice 
with Q = (independent ID chains) is also shown. Statistical 
errors are much smaller than the symbols. 

Here quantum Monte Carlo results for the cases n = 
1, 2, 3, and 4 will be presented. In addition, results will be 
shown for the conventional single-chain effective Hamil- 
tonian (^, which corresponds to ti = 0. This ID Hamil- 
tonian has been studied by Schulz via a mapping to a 
solvable fermion model in the continuum limit , with 
the result M ~ O.TlOy^ for the 2D system. The mapping 
has not been demonstrated rigorously, however, and the 
above form of M can furthermore be valid only for small 
a. Numerical calculations have previously been carried 
out for L < 10 [0, which is not sufficient for addressing 
the behavior for a <C 1 in the thermodynamic limit. 

Before presenting the mean-field results, quantum 
Monte Carlo calculations for the full 2D Hamiltonian 
(^ will be discussed. For the spatially isotropic system 
(a = 1), very accurate results for M have previously |2^] 
been obtained using ground state results for the staggered 
structure factor, 

N N 

i=i j=i 

Accounting for rotational averaging, the sublattice mag- 
netization is given by 

= SS{n, tt)/N {N -> oo). (6) 

Here this quantity will be extrapolated to infinite size for 
a < 1. Using the stochastic series expansion method with 
an efficient cluster update |2l[ , ground state properties of 
systems with several thousand spins can be studied (for 
the largest lattices considered here, inverse temperatures 
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P = J/T as high as 2048 were used in order to obtain 
results free of temperature effects) . 

For a = 1, the leading finite-size corrections to as 
defined in Eq. are positive and ~ 1/y/N This 
can be expected also for < a < 1 if the system is 
ordered. Fig. |l| shows results for a = 0.05 on L x L 
lattices with L up to 40. The results extrapolate to M > 
0, but subleading corrections to the linear behavior are 
clearly large. Previously, results for smaller L were used 
as evidence that M vanishes below a critical value ac ^ 
0.1 — 0.2 [^^. Results for rectangular lattices with 
different aspect ratios R — L^/Ly reveal a considerable 
dependence on R, as also shown in Fig. |l|. For R = 8, 
the expected linear behavior can be seen clearly, and for 
i? = 4 there is a cross-over to this behavior for large 
systems. For R = 2 there is a clear minimum, and the 
R = 1 results also suggest one. In the two latter cases 
the finite-size behavior is hence non-monotonic, and there 
has to be a maximum for even larger systems before the 
asymptotic, linear (with positive slope) approach to the 
infinite-size value, which for a = 0.05 is S{TT,Tr)/N 
0.0056 (from an extrapolation of the R = 8 data). 

The non-monotonicity can be understood as resulting 
from a cross-over from ID to 2D behavior. A chain 
of length Lx has an excitation gap A(La;) ~ i/L^- 
If this gap is larger than the effective energy scale of 
the coupling of the chains, i.e., the spin-stiffness p^, 
then the system essentially behaves as a system of ID 
chains, with exponentially damped correlations between 
the chains. A cross-over to 2D behavior can be expected 
when A (La;) ^ which occurs for smaller system sizes 
N = LxLy when the aspect ratio R is large, in agreement 
with the results in Fig. |^. When a ^ 1 (and therefore 
p| <C 1), quadratic lattices therefore have to be very 
large for extrapolations to infinite size to be meaningful. 
Instead, rectangular lattices with R increasing with de- 
creasing a should be used. Using aspect ratios as large as 
R = 16, the sublattice magnetization was calculated for 
a as small as 0.02. Below, the results will be compared 
with the single- and multi-chain mean field theories. 

In the single-chain theory the magnetization curve 
M{J±) is directly obtained from a calculation of M{h) 
for the Hamiltonian (||) . The effective model ^ depends 
explicitly on J±, however, and for each J± a search for the 
self-consistent values ft., Ai, . . . , A„ is required. With the 
Monte Carlo method used the derivatives dMj/dh 
and dMj/dXk can also be calculated. Using these, an 
iterative scheme where 

h{m+l) = h{m) + Ah{m), 
Afe(m + l) = Afe(TO) + AA,(m), (7) 

can be employed, starting from estimated values h{0) and 
Afc(O). The self-consistency conditions (^ give the correc- 
tions A/i(m) and A\^, (m) as the solution of n + 1 coupled 
equations, e.g., for n = 1, 
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FIG. 2. Self-consistent staggered magnetization vs inter- 
chain coupling in the single-chain mean-field theory (n = 0) 
and multi-chain mean-field theories with n = 1 — 4. Sta- 
tistical errors are at most comparable to the symbol sizes. 
Monte Carlo results for the full 2D Hamiltonian are shown 
with estimated error bars. The dashed line is the analyt- 
ical n = result [11]. The solid curves are of the form 
M„/Va = An(l + ba)ln'{a/a), with b = 0.095, a = 1.3, 
and 7 = 1/3 in all cases. These parameters, and Aq = 0.529, 
were chosen to fit the n = data. Only the amplitudes An 
were subsequently adjusted to fit the other data sets. 

Ah{dMo/dh-dMi/dh) 

+Ax, {dMo/dXi - dMi/dXi) = Ah ~ Mq, (8) 
Ah{dMo/dh + dMi/dh - 2/ Jj_) 

+Ax,idMo/dXi + dhh/dXi) = 2ft/ - Mi - Mq. 

Self-consistency is typically achieved this way in as few 
as two or three iterations. 

For a finite system, the self-consistent M vanishes be- 
low a critical value ac{L) which decreases with increas- 
ing L. In order to study the behavior for small a very 
large L have to be used. The largest sizes used here were 
L = 1024 for n = 0, 512 for n = 1, 2, and 256 for n = 3, 4. 
Inverse temperatures (3 = J /T as high as 2L were used 
in order to completely project out the ground state. 

All results for A/, including those for the original 2D 
Hamiltonian ([l]) extrapolated to infinite size, are shown 
divided by ^/a in Fig. |^. The behavior predicted by 
Schulz |ll| using a mapping of the n = mean-field 
theory to a solvable continuum model should then be 
a constant. The numerical results for n = do not 
agree with this; instead MqI-Jol appears to diverge as 
a ^ 0. The behavior for a < 0.4 is closely reproduced 
by the form A/q — ^o\/a(l -I- 6a)ln'''(a/a), with 7 « 1/3, 
Aq k, 0.53, a ~ 1.3, and h « 0.1. This result shows 
that the mapping of Eq. (|2|) to the continuum model is 
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FIG. 3. Self-consistent anisotropy parameters vs interchain 
coupling for n — 1,2, 3. All Xk, k — 1, . . . ,n, for given n are 
shown using the same symbols, and in all cases Ai < . . . < A„. 



not exact. A reason for this could be the presence of 
marginally irrelevant operators, which are known to lead 
to logarithmic corrections to physical observables in the 
case h = . The results for higher n also show a sim- 
ilar divergent behavior, but with the available computer 
resources it was not possible to extend the calculations 
to as small a as for n = 0. The above logarithmic form 
fits quite well also all the multi-chain results, with only 
the over-all factors An adjusted. This is a strong indi- 
cation that the logarithmic correction survives in the 2D 
limit (n — s- oo). The curves indeed approach the results 
obtained using finite-size extrapolations for rectangular 
2D lattices, confirming that the multi-chain mean-field 
theory converges correctly. Remarkably, the same ex- 
pression that describes all the mean-field data also fits 
the 2D results, with the amplitude A2D ~ 0.39. 

The self-consistent values of the xy-anisotropy param- 
eters are graphed in Fig. ^ for n = 1, 2, 3. For n> 1, the 
anisotropy is always largest at the edges, as expected, 
and rapidly decreases as the center chain is approached. 
The behavior for a ^ suggests a very slow asymptotic 
decay to zero — again an indication of log-corrections. 

To conclude, both the multi-chain mean-field theory 
and calculations for the original Hamiltonian strongly 
support a critical coupling ac = 0, and a staggered mag- 
netization that for small interchain couplings behaves as 
M ~ ^/a enhanced by a logarithmic correction. In the 
conventional single-chain mean-field theory {n ~ 0), all 
interchain quantum fluctuations are neglected. 2D quan- 
tum fluctuations develop systematically in the multi- 
chain theory as n is increased. For a ^ 1 , the functional 
form of the sublattice magnetization is the same for all 
n considered (n = — 4), indicating that the interchain 
quantum fluctuations only affect the over-all magnitude 
of M. Hence even the conventional single-chain theory 



gives the correct functional form for M, although the 
magnitude is over-estimated by a factor « 1.35. The 
previous analytical treatment of the single-chain theory 
pl| misses the log-correction. 

For the model considered here, it was possible to ex- 
plicitly test the multi-chain mean-field theory against 
large-scale quantum Monte Carlo results. In general this 
would not be possible, e.g., for systems with frustrated 
interactions where Monte Carlo simulations suffer from 
sign problems. The density matrix renormalization group 
method could be used to study the effective multi- 
chain models in such cases. 
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